/*=============================================================================
*
*   Copyright (C) 2020 All rights reserved.
*
*   Filename:		incremental_load.hpp
*
*   Author: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Date: 2020-03-05 17:17
*
*   Last Editors: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Last modified:	2020-08-04 01:43
*
*   Description:
*
=============================================================================*/
#ifndef INCREMENTAL_LOAD_H
#define INCREMENTAL_LOAD_H

#include "geometry/shape/elemental_dof.hpp"
#include "geometry/shape/shape_interface.hpp"
#include "model/solve/update_time_step.hpp"
#include <typeinfo>

namespace geoxfem {
namespace model {

struct incremental_load {};

template <>
bool check_convergence<incremental_load>::/*{{{*/
  apply(const VectorXR& ans, _real& tol) {
  _real norm     = ans.norm();
  _real relative = relative_tolerance::set_get();
  _real absolute = absolute_tolerance::set_get();
  /* if relative tolerance is set */
  _real tol1 = tol + relative * norm;
  /* if absolute tolerance is set */
  _real tol2 = absolute;
  /* take the larger one */
  tol = tol1 > tol2 ? tol1 : tol2;

  GEOXFEM_LOG_SEV(info) << "Tolerance: " << tol << ". Increment: " << norm;

  return tol >= norm;
} /*}}}*/

template <>
template <typename DimensionType, typename degree>
void update_tensor<incremental_load>::/*{{{*/
  apply(geometry::ShapeInterface<DimensionType, degree>& shp) {
  /* 1. update gauss points' status */
  VectorXR local_displacment = geometry::get_value_list_at_elemental_dof<
    displacement, dof_value_at_curr_iter>::apply(shp);
  const material::MaterialInterface<DimensionType>* mat = shp.getMaterial();
  for(auto& gp : gpsIn(shp)) {
    /* calculate epsilon increment */
    VectorXR voigt = BeeMatrix(shp, gp) * local_displacment;
    arithmetic::Strain<DimensionType> dde;
    dde.fromVoigt(voigt);
    gp.addEpsilon(dde);
    /* calculate sigma increment */
    arithmetic::Stress<DimensionType> old_sigma = gp.getCurrSigma();
    arithmetic::Stress<DimensionType> new_sigma = old_sigma;
    util::iRecord<DimensionType>&     param     = gp.getParam();
    mat->IntegralThroughEpsilonPath(dde, new_sigma, param);
    arithmetic::Stress<DimensionType> delta_sigma = new_sigma - old_sigma;
    gp.addSigmaAtTimeStep(delta_sigma);
  }
  /* 2. mapping new gauss points' status to vertex*/
  algorithm::G2V_mapping(shp);
  return;
} /*}}}*/

template <typename CouplePolicy>
struct update_iteration<CouplePolicy, incremental_load> { /*{{{*/
  template <typename Mesh, typename DOFHandler, typename TimeStep>
  static void
  apply(const Mesh& mesh, const DOFHandler& dof_handler, TimeStep& tstep) {
    GEOXFEM_LOG_SEV(info) << "Updating within time step";
    size_t topo_dim = traits::get_topology_dimension<Mesh>::apply(mesh);
    /* 1. update tolerance at current iteration */
    tstep.m_converged =
      check_convergence<incremental_load>::apply(tstep.m_answer, tstep.m_tol);
    /* 2. update unknowns at current iteration */
    for(auto dof : dofsIn(dof_handler)) {
      _real value = tstep.m_answer[dof->getIndex()];
      dof->setValueAtIteration(value);
      dof->addValueAtTimeStep(value);
    }
    /* 3. reset gauss on vertices */
    for(auto& vertex : shapesIn(mesh, 0))
      vertex->toVertex().getGauss().reset();
    /* 4. update stress/strain at gauss */
    if(traits::if_consider_deform<CouplePolicy>::value)
      for(auto shp : shapesIn(mesh, topo_dim))
        update_tensor<incremental_load>::apply(*shp);
    return;
  }
}; /*}}}*/

template <typename CouplePolicy>
struct update_time_step<CouplePolicy, incremental_load> { /*{{{*/
  template <typename Mesh>
  static void apply(const Mesh& mesh) {
    GEOXFEM_LOG_SEV(info) << "Updating between time step";

    /* 1. update unknowns at current time step */
    for(auto vertex : shapesIn(mesh, 0))
      for(auto& dof : dofsIn(vertex->toVertex())) {
        dof.second.addValue(dof.second.getValueAtTimeStep());
        dof.second.setValueAtTimeStep(0);
      }
    /* 4. update stress/strain at gauss */
    size_t topo_dim = traits::get_topology_dimension<Mesh>::apply(mesh);
    for(auto shp : shapesIn(mesh, topo_dim))
      for(auto& gp : gpsIn(*shp)) {
        gp.addSigma(gp.getSigmaAtTimeStep());
        gp.setSigmaAtTimeStep(Matrix3R::Zero());
      }
    return;
  }
}; /*}}}*/

}  // namespace model
}  // namespace geoxfem

#endif /* INCREMENTAL_LOAD_H */
